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O ' Abstract 

<D ■ 

. In this paper we present a method for estimating unknown parameters that appear on an avascular, spheric 

tumour growth model. The model for the tumour is based on nutrient driven growth of a continuum of live cells, 
whose birth and death generate volume changes described by a velocity field. 

The model consists on a coupled system of partial differential equations whose spatial domain is the tumour, 
that changes in size over time. Thus, the situation can be formulated as a free boundary problem. 

After solving the forward problem properly, we use the model for the estimation of parameters by fitting 
the numerical solution with real data, obtained via in vitro experiments and medical imaging. We define an 
appropriate functional to compare both the real data and the numerical solution. We use the adjoint method 
for the minimization of this functional, getting a better performance than the obtained with the pattern search 
method. 

Key words: avascular tumour, PDE constrained optimization, inverse problem, mathematical modeling, adjoint 
method 

1 Introduction. 

The interest for research in modeling cancer has grown enormously over the last decades (2j[rj|5j|6l, and it became 
one of the most challenging topics involving applied mathematicians working with researchers in the biological 
sciences. One of the main motivations is the fact that, according to the World Health Organization, about six million 
people die annually because of cancer, being this one the second main fatal disease in the industrial countries. 

Key comments on the importance of mathematical modeling in cancer can be found in a vast part of the liter- 
ature. For example, in the work by Bellomo et al. (61, they emphasize the fact that "applied mathematics may be 
able to provide a framework in which experimental results can be interpreted, and a quantitative analysis of external 
actions to control neoplastic growth can be developed". Moreover, "models and simulations can reduce the amount 
of experimentation necessary for drug and therapy development". 

In this paper we consider the case of avascular multicellular spheroids (MCS). Pioneers in this subject have 
been, for example |[T2ll20ll . where the first spatio-temporal models of MCS' growth have been developed. The study 
of MCS is interesting because they provide the best insight into the effects of varying nutrient concentrations or 
the effectiveness of chemotherapeutic drugs on tumours in vivo, and their behaviour can be studied experimentally 
(in vitro) by controlling environmental conditions in which they grow: for example, the radii of the tumour can be 
monitored while changing the chemotherapeutic drug or oxygen levels. 

In addition, another variables can be measured. If possible, experimentalists can get information about the 
distribution of substances within the tumour. Moreover, via medical imaging, histopathology and potentially other 
sources, they can also get data about the density of the different kind of cells conforming it: proliferating, quies- 
cent, necrotic. For instance, as documented in |fl9l the Boron Neutron Capture Therapy (BNCT) technique gives 
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information about the evolution in size of a melanoma, and in lfl4l . they obtain information about the growth of a 
glioma via Magnetic Resonance Imaging (MRI). 

That is why in this general approach of modeling, the key variables are the tumour size (radius), the concen- 
tration within the tumour of growth-rate limiting diffusible chemicals (nutrients such as oxygen or glucose or a 
chemotherapeutic drug) and the density of cells. Since the tumour changes in size over time, the domain on which 
the models are formulated must be determined as a part of the solution process, giving a vast class of moving 
boundary problems ll9l[T0l. 

In this article, we propose a framework for estimating unknown parameters via a PDE-constrained optimization 
problem, following the PDE-based model by Ward and King li23l . and considered also by Knopoff et al. ifToll . In 
this approach, avascular tumour growth is modeled via a coupled nonlinear system of partial differential equations, 
which makes the numerical solution procedure quite challenging. 

This kind of problem constitutes a particular application of the so-called inverse problems, which are being 
increasingly used in a broad number of fields in applied sciences. For instance, problems refered to structured 
population dynamics lfT8l . computerized tomography and image reconstruction in medical imaging |[2T1 l24l . and 
more specifically tumour growth (3l[l4l[T6l, among many others. 

Extending the work done in lfT6l . we are concerned with developing a robust PDE-constrained formulation that 
let us find the best set of parameters of a tumour growth model that fits patient or experimental data. In contrast 
to the previous work, in which we used a free-derivative method (pattern search algorithm) we now use the adjoint 
method in order to find the derivative of the afore mentioned functional. We want to find the parameters that would 
be of interest by defining a functional to be minimized. In this way, we would obtain the best set of parameters that 
fits patient-specific data. 

The contents of this paper, which is organized into 9 sections and an appendix, are as follows: Section 2 consists 
in some preliminaries about the model and the definition of the direct problem. It can be considered as a revision 
of our previous work. Section 3 deals with the motivation and formulation of the minimization problem. Section 4 
introduces the adjoint problem, deriving the optimality conditions for the problem. Specifically, we show how the 
adjoint method may be used to find the derivative of the solution of a PDE with respect to a parameter that does 
not appear explicitly in the equation. Section 5 refers to specification of the tools used in Section 4 for the concrete 
problem. Section 6 deals with the numerical solution of the adjoint problem, designing a suitable algorithm to 
solve it. In particular, we develop a method to deal with some singularities in the PDEs. Section 7 deals with the 
minimization method to be used and, specifically, we propose an algorithm to find a minimum of the functional, 
by choosing a proper descent direction using derivatives. In section 8 we show some numerical simulations to 
give information on the behaviour of the function and its dependence on the parameters. Section 9 presents the 
conclusions and introduces some future work related to the contents of this paper. 

Some words about our notation. We use (•, •) to denote the L 2 inner product (the space is always clear from the 
context) and we consider the sum of inner products for a cartesian product of spaces. For a function F : f x W — s> Z> 
such that ((j),p) F(§,p), we denote by F'(§,p) the full Frechet-derivative and by §x(<|>,/>) and |^((j),p) the partial 
Frechet-derivatives of F at (§,p). For a linear operator T : t y — > Zwe denote T* : Z>* — > 'Jf* the adjoint operator of 
T. If T is invertible, we call T~* the inverse of the adjoint operator T*. 

2 Some preliminaries about the model. 

In lfl6l we dealt with the resolution of a coupled system of spatio-temporal PDEs which involves initial and bound- 
ary conditions, with the additional difficulty that the boundary is also an unknown. In that work we considered 
that the tumour is a spheroid which consists of a continuum of living cells, in one of two states: live or dead. The 
rates of birth and death depend on the nutrient and chemotherapeutic drug concentration. It is supposed that those 
processes generate volume changes, leading to cell movement described by a velocity field. In order to develop 
a method capable to recover parameters, we will consider the simpler case in which no treatment is supplied to 
the neoplastic formation, leaving this case for future work. Under this assumption, the system of equations to be 
studied is: 
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l + i^l = M5.e)-M«.e)h. C") 

1 3(r 2 v) 



- [V L * m (q,e)-(V L -V D )* rf (q,e)]Ti, (2.3) 

r z or 

where the dependent variables r|, 5 and v are the live cell density (cells/unit volume), nutrient concentration and 
velocity, respectively The independent variables are the radial position r inside the tumour and time t. Constants Vl 
and Vd correspond to the volume of a living and a death cell, respectively. The number D is the diffusion coefficient 
of the nutrient and (3 is a positive constant related to the nutrient's consumption rate. As it is described in |23l . 
equation (12.1b states that the rate of change of r\ is dependent on the difference between the birth k m (q, 6) and death 
k c /(q, 6) rates (6 is a vector of parameters associated to these functions). The functions k m and kd are taken to be 
generalized Michaelis-Menten kinetics with exponent 1, i.e. 



fc m M) = A I —-j— I , (2.4) 

ka(q,B) = fi(l-o-i-Y (2.5) 

where 6 = [A,B,<; c ,<^/,a] r are model parameters. As stated in the appendix of |[22l . there appears to be no appro- 
priate data available on the parameters C4 and a, constituting one extra motivation for this work. 

Inherent in this problem are two timescales: the tumour growth timescale (?» 1 day) and the much shorter 
nutrient diffusion (ps 1 min), letting us to adopt a quasisteady assumption in the nutrient equation (see |[22l ). 
Therefore, we replace (12.21) by the quasisteady approximation 

1 d ( 2 dq\ p 



2.1 Initial and boundary conditions. 



r l f)=^k m (q,Q)T\. (2.6) 



As it has been mentioned, the tumour is assumed to be a spheroid that exhibits radial symmetry. That is why, not 
only the state variables r|, q and v are important, but also the tumour radius is a key variable to be determined. 
Since the tumour changes in size over time, the domain on which the model is formulated (and the PDEs are valid) 
must be determined as part of the solution. 

Let 5(0 be the tumour radius at time t. At t = we will consider the tumour at a certain stage of its evolution. 
Hence the initial conditions are a known radius 5(0) and an initial live cell density 

n(r,0) =Tu(r). 

Because symmetry is assumed about the tumour center, there is no flux there. That is why the boundary conditions 
at r = are: 

|r (0,0 = 0, (2.7) 

or 

v(0,/) = 0. (2.8) 

Moreover, on the external boundary (which is also the boundary of the complement of the tumour as a subset 
of the body), the following conditions are taken: 

S (5(0,0 = co, (2.9) 

d 4; = v(5(0,0, (2-10) 
dt 

where cq is the external nutrient concentration. 
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2.2 Nondimensionalization and fixed domain method. 



Following the ideas exposed in El [T] |22l [23l, the mathematical model is rescaled and the domain [0,5(0] °f 
the tumour is transformed onto the interval [0, 1]. This is a very useful approach when dealing with free-boundary 
problems, as mentioned in flOl . Hence, let us define the following functions 

N(y,t) = V L r\(yS(t/A),t/A), 
C(y,t) = -q(yS(t/A),t/A), 

V(y,t) = J- V (yS(t/A),t/A), 
Ar 

S(t) = -S(t/A), 
a(c,*) = j[k m (c,^)-k d (c,^)}, 

b{c,V) = j[k m (c,-&)-(l-8)k d (c,-&)], 
*(c,0) = j& m (c,fl), 

where r = (3V l /4k) 1 / 3 is the radius of a single cell, 8 = V D /V L , P = r$/(V L coD) and ■& = \A,B,c c ,ca,a] with 
c c = Cx/cq and c d = C^/cq. Thus, we obtain the following system to be solved: 



Nt-~yNy + ^N y = N[a(C,V)-b{C,$)N], 0<y<l, t>0, 
2 



Cyy + -Cy = ^C^S^N, < J < 1 , ?>0, (2.12) 

2 

V y + -V = b(C,$)NS, 0<y<l, t>0. (2.13) 



The initial conditions for the transformed problem are: 



N(y,0) = N r (y), 0<y<l, (2.14) 
S(0) = S 7 , (2.15) 



where Nj(y) = VifliiyS (0),0) and Si = 5(0)/ro, and the boundary conditions are: 



V(0,f) = 0, t>0, (2.16) 

C y (0,t) = 0, t>0, (2.17) 

C(l,f) = 1, t>0, (2.18) 

S'(t) = V(l,t), t>0. (2.19) 



From now on, equations (12.11l )-( |2~T9l ) will be referred to as the direct problem. 



3 Formulation of the minimization problem. 

As described above, there is a set of parameters (some of them unknown) that determines the behaviour of a tumour 
growth. For this reason we propose to use an inverse problem technique in order to estimate them. 
We define the following vectors: 

4 = [N,V,C,S] T , (3.1) 
P = [c c ,Cd,G} T , (3.2) 
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Figure 1 : Microscopic image of neoplastic colonies that grow with an external nutrient supply. Courtesy of CNEA 
( Comision National de Energi'a Atomica). 



where (]) represents the solution of the direct problem (the components of (j) are the state variables of the problem) 
for each choice of the vector of parameters r> = (A,B,p), where A and B are assumed to be constants. Hence from 
now on, we will use just p instead of t> as the vector of parameters. 

Let us assume that experimental information is available during the time interval < t < T. Then, the general 
problem we are interested to solve can be formulated as: 

Find a parameter p able to generate data § = [N,C,V,S] T that best match the available (experimental) 

information over time <t <T. 

For this purpose, we should construct an objective functional which gives us a notion of distance between the 
experimental (real) data and the solution of the system of PDEs for each choice of parameters p. 

First of all, it is important to decide which variables are capable to be measured experimentally. For instance, it 
is clear that the tumour radius can be known at certain times tk,k= 1, ...,M via MRI, PET (Positron Emission To- 
mography) or CT (Computed Tomography). For example, Figure [His a microscopic field that shows the formation 
in vitro of neoplastic colonies which grow as spheroids with an external nutrient supply. Such experiments could 
help to determine optimal variables and parameters in order to control real tumour growth. 

So, the first possibility for defining a functional could be: 



where S(t) is the radius evolution obtained by solving the direct problem for a certain choice of p and S*(t) = 
jrS*(t/A) is the evolution measured experimentally (real data). 

Another variable that could be measured is the density of living cells, also via biomedical imaging. As it 
was mentioned before, this could be done via PET technique for a tumour in vivo, or via immunofluorescence 
and electronic scan microscopy technique for in vitro cases. Thus, we are motivated to define a functional that 
reproduces in a better way the knowledge we have about the process: 




(3.3) 



J(N,S,p) = ^ f 1 f T [N(y,t)-N*(y,t)] 2 dtdy + ^ C [S(t) - S*(t)] 2 dt 
2 JO JO 2. Jo 



(3.4) 
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where N(y,t) and N*(y,t) are the living cell concentrations for the direct problem solved with the parameters p 
and the real data, respectively (both of them in the domain [0, 1] x [0, 7 1 ]). The positive constants /ji and ^2 are 
introduced, as we shall see, to take into account the different order of magnitude between N and 5. Note that, for 
instance, if we take ^1=0 and ^2 = 1 in ( 13.41 ) we get ( 13.31 ). In this way, these two parameters will give us some 
flexibility in order to choose an appropriate functional according to the experimental method used to obtain the 
data. 

It should be noted that the spatial integration is done over the interval [0, 1], because we are using the solution 
in the fixed domain. 
Let us define 

N -N y ^-y + ^N y - N (a(C, p) - b(C, p)N) 



2 

Vy + -V 

' y 



c 



YV 



2 -c, 

y ' 



■b(C,p)NS 
-k(C,p)NS 2 



(3.5) 



V(l,-)-S> 

no,-) 

C(l,-)-l 

Cy(0,-) 
N(;0)-N! 

5(0) -5/ 

In this way we can rewrite the system of PDEs (12. 1 lb - d27T9b described in the previous section as E(§,p) = 0. 

The set of parameters that best matches the experimental data with the generated data provided by the direct 
problem can be computed solving a PDE constrained optimization problem, namely: 



minimize J($,p) 
p 

subjectto E($,p) = 0, 
P G Uad, 



(3.6) 



where U a d denotes the set of admissible values of p. In our case, according to ( 13.21 ), U a d should be a subset of R 3 . 
Notice that a solution (§,p) must satisfy the constraints E(§,p) = 0, which constitute the direct problem. 

We remark that, in general, there is a fundamental difference between the direct and the inverse problems. 
In fact, the latter is usually ill-posed in the sense of existence, uniqueness and stability of the solution. This 
inconvenient is often treated by using some regularization techniques ll2Tl[TTl[T5l . 



4 Formulation of the reduced and adjoint problems. 

In the following, we will consider a generic optimization problem, which has the form: 

minimize J($,p) 
p 

subjectto E(ty,p) = 0, C 4 - 1 ) 

P € U ad , 

where / : y X U a d — > R is an objective function and E : y x U a d — > Z is a state equation, for y and Z Banach 
spaces and U a d is a set of admissible points. 

For completeness, in this section we will present a general theory in order to solve problem (14- 1 b - According to 
the ideas exposed in (HEl, we make the following assumptions: 

(Al) U a d € R" 1 is a nonempty, closed and convex set. 
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(A2) J : y x U a d M and £ : 9^ x £/ c „i — > Z are continuously Frechet-differentiable functions. 

(A3) For each p G ?7^/ there exists a unique corresponding solution (()(/?) G J" such that £((()(/;>),/?) = 0. Thus, 
there is a unique solution operator p G U a d i— > (j)(p) G 9^- 

(A4) The derivative : 9^ — s- .Z is a continuous linear operator, and it is continuously invertible for all 

P G U cl d- 

Under these hypotheses §(p) is continuously differentiable on p G U a d by the implicit function theorem. Thus, 
it is reasonable to define the following so-called reduced problem 

minimize J(p)= J(§(p),p) 

p (4.2) 
subject to p G U a d, 

where (j)(p) is given as the solution of E{§(p),p) = 0. 

In order to find a minimum of the continuosly differentiable function J, it will be important to compute the 
derivative of this reduced objective function. Hence, we will show a procedure to obtain f by using the adjoint 
approach. Since 

(J'(p),q) = (^^{p),pU'(p) q ^ + (^(^p),p),q 
^(p)Y^p),p) + ^(p),p),q 

f'ip) = ^(p)Y^(p),p) + ^{p),p). (4.3) 
Let us consider X G Z* as the solution of the so-called adjoint problem: 

^(p),P)+{^(m,P)J^ = 0. (4.4) 

where (jj§(§,p)j is the adjoint operator of ||((]),/?). Note that each term in ( 14.41 ) is an element of the space 9'*- 
An equation for the derivative §'(p) is obtained by differentiating the equation E{§{p),p) = with respect to 

p: 

dE dE 

— (^p),pW(p) + —(^p),p)=0, (4.5) 

where is the zero vector in Z. 
By using (14.31 ) we have that: 

j'(p) = ^Xp)y d ^p),p) + d J-^p),p) 

d ^(p), P )\\+^(p),p), 



we see that 



where in the second equation we used (14.51 ) and for the last equation we used (14.4I ). Then: 

J'(p) = ^{p),p)+ (^(p),p)Yx. (4.6) 

Notice that in order to obtain J'(p) we need first to compute §(p) by solving the direct problem, followed by 
the calculation of X by solving the adjoint problem. For computing the second term of (14.61 ) it is not necessary to 
obtain the adjoint of ^(<^(p),p) but just its action over X (see the Appendix). 
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5 Getting the adjoint equation for the concrete problem. 



For our case, let us define CI = [0, 1] x [0, T], the spatio-temporal domain of interest. Note that, according to (13.ll ). 
(j) is an element of a suitable vector space. Let us consider the function spaces 

y={c l (Q.)) 2 xC 2 (n)xc l ([o,T}), 

Z = (C 1 ^)) 2 x C 2 {Q.) x {C l ([0,T])) 4 x C([0, 1]) x R. 

The spaces C l and C 2 inherit the inner product from L 2 , so the spaces f and Z are Hilbert spaces (therefore we 
can identify < J r * an d Z>* with y and Z respectively). It is worth mentioning that we consider these vector spaces 
because we look for strong solutions of the PDEs, i.e., we require differentiability of the state variables. 

In order to obtain the adjoint operator of we have to find f |? J such that: 



dE \\ 



(5.1) 



Hence, we define the directions n, v, c, and s for the state variables N, V, C and S. Let g = [n, v,c,s] T , then 



dE 
d4» 



■ lim 



E($ + yg,p)-E(§,p) 



After some algebraics, it can be shown that M (§,p)g is given by: 



«? + ^—i—ny 



s'S-S's ; 



^N y y+N y *^-[a-bN\n-N %c-%Nc-bn 



da 



db 



Vy + -V ■ 



3fo 
3C 



c y.v + y c y 



NSc - bSn - bNs 



(5.2) 



■ kS 2 n - ^NS 2 c - 2kNSs 

v(l,-W 
v(0,-) 

c,(0,-) 
n(-,0) 
,(0) 

Note that E(<ty,p) and X should have the same number of components. Also, each component of X must be in 
a subspace of the corresponding component of E. For example, the first three components of X must depend on 
space and time, the fourth one only on time, the last one is just a real number, and so on. 

So we define: 



X( y ,t) = [X 1 ( y j)M(y,t),h{y,t)M(t),h(t)M(t),h(t)M(y)M] T - 



(5.3) 



An inspection over equations d5.ll ) and (15.21 ) shows that, roughly speaking, we should remove the spatial and 
temporal derivatives from g and pass them to X. 

The calculations make use of successive integration by parts to express each derivative of g in terms of a 
derivative of X. Omitting here the details, that are shown in the Appendix, we obtain the following system of 
equations, which constitutes the adjoint problem (14.41 ): 
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-h, - (^y^) ~ +a ~ 2bN ) Xi ~ bSkl ~ kSlh =* i i(N*-N), (5.4) 

X 2y --X 2 -^X 1= 0, (5.5) 

2, (2 dk ,\„ 3& „ (da db \ „ 

X 3> , - + ^ - j X 3 - ^ -iV (- - X, = 0, (5.6) 

X l (y,T)=0, (5.7) 

^(1,0 = -U(t), (5.8) 

X 3> ,(0,0=0, (5.9) 

^(1,0=0, (5.10) 

Ut(t) = f Q (^h-~{Nyh) + bNX 2 + 2kNSh^dy + f i 2 (S*(t)-S(t)), (5.11) 

MT) = -j Q — ydy, (5.12) 

^(0=^(0,0, (5-13) 

X 6 (0=X 3y (l,0 -2X3(1,0, (5-14) 

X 7 (0 = 3X3(0,?), (5.15) 

X 8 (y)=X 1 (y,0), (5.16) 



5(0) 



rf;y-X 4 (0). (5.17) 



Equations (15.4I )- (15.17I ) shall be solved in order to get X. Notice that the adjoint equations are posed backwards 
in time, with a terminal condition at t = T, while the state equations are posed forward in time, with an initial 
condition at t = 0. 

It is worth stressing that an explicit expression for X4 can be obtained. In fact, taking into account that Xi (y, T) = 
0, we have that 

Then, by equation [2. 1 1 I from the direct problem, we have that 

N yjy-^ N y = N t- N i a - bN }- 
But using the fact that X 4 (T) = 0, then — X 4 (0 = f t T X^dx, i.e., 

r T r l (y d NvV \ r T 

X 4 (0 = J J o i^j t (Nyki ) - -^-Xi - bNX 2 - 2kNSX 3 J dydx + fi 2 J (S-S*)dx 

= f £{[N t -N(a- bN)} ^ - bNX 2 - 2kNSl 3 } dydx - £ y Ny(y '*^ M dy + ix 2 ^ (S - S* )dx. 



6 Designing an algorithm to solve the adjoint problem. 

It is worth stressing that obtaining model parameters via minimization of the objective functional / is in general 
an iterative process requiring the value of the derivative. To compute /' we just solve two systems of PDEs per 
iteration: the direct and the adjoint problems. This method is much cheaper than the sensitivity approach lfT3l 



9 




Adimensional time 



Figure 2: Evolution of the tumour radius in time. 

in which the direct problem is solved many times per iteration. That is why we developed an implementation in 
Fortran 2003 using an object-oriented strategy (with Fortran Intel Compiler 12.0.3). For the direct problem, Figure 
|2]shows the evolution of the tumour radius in time, and Figure [3]represents the living cell density within the tumour 
for two different times. 

Although the adjoint problem is quite similar to the direct one, there are more difficulties to solve it. For 
example, there is no explicit boundary condition for %2- In our particular case, it is not necessary to compute {X,;} (=4 
in order to calculate the derivative of / with respect to p, because {Ej} i=4 does not depend on the parameters p (see 
equations (13.51 ) and (14.6I )). However, X4 is required since it gives us the boundary condition for X2 (see equation 
( 15.81 )). To design a numerical procedure we perform the following steps at time T: 

- Equation ( 15771 ) states that ki (-,T) = 0. 

- By equation (15.121 ) we have that ^(r) = 0, which gives us a boundary condition for X2 (see equation (15.81) ). 

- Equation (15.51 ) can be solved analytically getting A^(-, T) = 0. 

- Equations (15.61) . ( 15791 ) and (15.101) allows us to obtain X 3 (-, T). 

Knowing the solution at time t, we obtain the solution at time t — At in the following way: 

- By equation ( 15.41 ) we first obtain X\ t (•,/). Then we get 'k\{-,t — At) using a backward finite difference. 

- Using equation ( 15.111 ), we integrate numerically to obtain X^t{t) and then we get \$(t — At) by means of 
backward finite differences. 

- With the value of X^{t — At) we obtain hi{\,t — At) via equation ( 15.81 ). 

- Equation ( 15.51 ) can be solved numerically to get 7^i{-,t — At). 

- Solving equations (15761) . d579l ) and (15.101) we obtain Xj(-,t — At). 



As well as in the direct problem, in the adjoint one we have to be careful with the singularities in the PDEs. For 
example, if we take a look to equation (15.6b together with the boundary conditions (15.91) and (15.101 ). for a fixed time 
t, we can ask ourselves about the solvability of this problem around y = 0. However, there is a difference between 
the direct and the adjoint problems regarding to the kind of singularities that equations (12.121) and (15.61 ) have in the 
origin. 
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Figure 3: Live-cell density within the tumour for two different times. Blue line corresponds to t = and red line 
corresponds to t = 2.5. Note that from t = to t = 2.5 the tumour has also grown in size. 



The second term in (15.6I ). for instance, looks harmless because Xs y (0,t) = by (I5.9I ), so upon expanding by 
Taylor about and dividing by y, the singularity disappears. On the other hand, the problem with the third term is 
harder, because we get a blowup in the origin. To solve this problem we transform equations ( I5.6I ), ( I5.9I ) and ( I5.10I ) 
into a first order ODE for a fixed time t, namely: 



v 



2 dk 



N' 



, db da 



dC 



= 0, 
v(0)=0, 



(6.1) 

(6.2) 
(6.3) 



where u(y) = X^(y,t) and v(y) = Xs y (y,t). Then, for a fixed £ > we propose a parameter q = v(l), and solve 
the system (16.1I )- (16.3I ) in the interval [e, 1] with boundary conditions u{\) = and v(l) = q, obtaining a solution 
[u q , v q ] T . Using Taylor expansions near y = we extend these solutions to the whole interval [0, 1] (see |4ll). 



u q (0) u q (e)-eu! q (e), 
v q (0) w v q (e)-ev' g (e). 



(6.4) 
(6.5) 



The next step is to define a function 

F(q)=v q (0), (6.6) 

and to find a root of F, i.e. to find q such that F(q) = 0. Then, the solution of the system will be [u q ,v q ] T extended 
to the interval [0,1]. 

To solve the first order ODE (15.5b for X2(-,t) with boundary condition A,2(l,f) known from (15 - 8b . we also 
solve the problem in the interval [e, 1] and then extend the solution to the interval [0, 1] using a first order Taylor 
expansion. 

In general, the derivatives that appear in the adjoint system of PDEs are approximated using a finite difference 
scheme. For example, in order to solve equation (15.41) we consider 



ki (y,t — At) m ki (y, t) - h t (y,t)At, 
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and using (15.41) we get 



M y,t-At) « W,t) + Ai(Wy + YMy, 

+ A* + +fl(C(y,0) -2fc(C(y,0 W,0) 

+ A? {fe(C(y, 0)5(0^2 (y, t) + k(C(y,t))S(t) 2 h (y, t) + N* (y,t) - N(y, t) } . 

7 Optimization. 

It is well-known ifTTl that gradient-based optimization algorithms require the evaluation of the gradient of the 
functional. One important advantage of evaluating the gradient through adjoints is that it requires to solve the 
adjoint problem only once per iteration, regardless the number of inversion variables. Note that the derivative of 
the functional can be approximated by using finite differences, but this is an expensive approach because it needs, 
for each optimization iteration, to solve the forward problem as many times as inversion variables are. 
The method we will use for minimizing the functional / can be summarized as follows: 

Algorithm 7.1 Adjoint-based minimization method. 

1. Give an initial guess p° for the vector of parameters. 

2. Given the vector p k in step k, solve the direct and adjoint problems at this step. 

3. Obtain the derivative of the functional, i.e. J > (p k ), using \4.6\ . 

4. Move in the direction of —f (p k ), i.e., compute p k+l = Tlu ad \p k — oJ 1 (p k )\, where CC is a positive parameter 
to be chosen, and Tlu ad denotes the projection on the set of admissible points. 

5. Stop when J(p k+l ) is less than a tolerance TOL\ > 0, or when the distance between two consecutive itera- 
tions is less than a tolerance TOL2 > 0, that is, || p k+l — p k ||< TOL2 

8 Numerical experiments. 

The goal of this section is to test and evaluate the performance of an adjoint-based optimization method, by exe- 
cuting some numerical simulations of Algorithm [7j] for some test-cases. 

The living cell density and the tumour radius are generated via the forward problem. We show here the results 
obtained by assuming standard values c c = 0.1, c c i = 0.05, a = 0.9, as suggested in |[22l . 



8.1 Model-generated data. 

Consider first an optimization problem that consist in minimizing the functional (13 -4b . where N*(y,t) and S*(t) are 
generated via the forward model, for a choice of the model parameters c c = 0.1, c c i = 0.05, a = 0.9. 

The tumour is first detected at time t = 0, by which time it has grown following the model |[22l . Originally, at 
an adimensional time t — 4.5, a single cell started to take nutrients from the environment, letting it grow up to a 
dimensionless size Sj « 34. Thus, the initial profile for this test case is the one shown in figures [2] and [3] 

In order to work with functional (13 -4b . it is necessary to define the landmark points yj and the times tk where the 
measurements are made. For simplicity, and to be consistent with the way we solved the direct problem, we took 
the same spatial grid for the landmarks, i.e., 30 equidistant points = y\ < ... < 3^30 = 1. Regarding to the time 
selection, it is apparent from the experiments that taking the adimensional time T = 0.5 is sufficiently representative 
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Cd value 

Cc value 

Figure 4: Functional value J in terms of c c and Cd for constant a = 0.9. Note that the surface reaches a minimum 
near c c = 0.1 and c c i = 0.05 
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Figure 5 : Evolution of the functional value with the number of iterations. 

(it corresponds to 50 time steps of length 0.01). The factors p\ and ^ are taken to be 100 and 1 respectively, and 
the parameter a used in the projection over the admissible set is taken to be 0.1. 

Figure |4] shows the value that the functional (13.41 ) takes for different values of c c and c c i, remaining a as a 
constant. It is worth mentioning that J looks convex and that the variations are greater with respect to c c compared 
to those with respect to cj. 

The idea of this test case is to investigate how close the original value of the parameter can be retrieved. 
However, it is not a trivial one, because we do not know, for instance, if the optimization problem has a solution or, 
in that case, if it is unique or if the method converges to another local minima. 

We emphasize that we have run the algorithm several times using different initial random conditions and in all 
cases the results were similar. They can be summarized as: 

- Stopping criteria: functional value lower than 10~ 6 or norm of the gradient lower than 10~ 12 

- Iterations/elapsed time: 140/35 min 

- Initial point: p = [0.16,0.03,1.0] 

- Final point: p f = [0.1006492,0.084465653,0.9297853] 

- Functional final value: J(p f ) = 0.991496220 x 10~ 6 

Figure [5]represents the evolution in the value of / with the number of iterations. It is worth stressing that, even 
the stopping criteria required 140 iterations, taking just about 90 iterations would be sufficient to obtain similar 
results. In fact, Figures [6] and [7] show the evolution of c c and a respectively, and the real value of this parameters. 
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Figure 6: Evolution of the c c value with the number of iterations (blue line) and real value of c c (red line). 




Figure 7: Evolution of the a value with the number of iterations (blue line) and real value of a (red line). 



8.2 Model-generated data with random noise. 

It is well known that the presence of noise in the data may imply the appearance of strong numerical instabilities in 
the solution of an inverse problem Q. 

The outputs of the detectors and experimental equipments where the variables N* and S* are measured are often 
affected by perturbations, usually random ones. As stated in [3J, it is in general valid to consider a 5% of random 
noise. 

Once again, it is assumed that the tumour is first detected at time t = 0, by which time it has grown following the 
model ll22l . Originally, at an adimensional time t — 7, a single cell started to take nutrients from the environment, 
letting it grow up to a dimensionless size 5/ « 53. 

The landmark points, the adimensional time T = 0.5 and the factors /U2 and a are taken to be equal as in 
the previous case. Note that, although this requires more iterations and consequently more computational time, we 
would get more information from a tumour that was detected a little bit later. 

After running the algorithm several times using different initial random conditions, the obtained results were 
similar. They can be summarized as: 

- Stopping criteria: functional value lower than 10~ 6 or norm of the gradient lower than 10~ 12 

- Iterations/elapsed time: 140/35 min 

- Initial point: p = [0.08,0.07,0.93] 

- Final point: p f = [0.1105396,0.07723431,0.9172613] 

- Functional final value: J(p f ) = 0.01830000000 
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Figure 8: Functional value / in terms of c c and cj for constant a = 0.9, for data with a 5% of random noise. 




Figure 9: Evolution of the functional value with the number of iterations, for data with 5% of random noise. 



Figure [8] shows the value that the functional (13.41 ) takes for different values of c c and cd, remaining a as a 
constant and assuming that 5* and N* are obtained with a 5% of random noise. 

Figure [9] represents the evolution in the value of J with the number of iterations. A comparison with Figure 
[5] shows that the functional values are greater in this case, but the algorithm stops because the variations become 
small. Figure [10] shows the evolution of c c and the real value of this parameter. 

We can choose one of the variables considered in the functional (13.41 ) and look for difference between the real 
value of this variable and the value that conesponds to the solution of the direct problem for the parameters obtained 
after running the algorithm. For example, let us take the tumour radius at time T. Figure [TT] shows a sequence of 
spheroids obtained in some of the iterations. As the difference between the radii is small compared to the radii 
themselves, we make a zoom, obtaining Figure [ 



8.3 Comparison between the two cases with and without noise. 

The case in which we considered model-generated data with 5% of random noise is, as expected, not as precise as 
the case in which the data is generated without noise. 
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Figure 10: Evolution of the c c value with the number of iterations (blue line) for data with 5% of random noise and 
real value of c c (red line). 




Figure 11: A sequence of spheroids obtained by Algorithm [7j] and choice of one point in the boundary of the real 
boundary and other point in the boundary of the tumour obtained after 60 iterations. 




Figure 12: A zoom from Figure QT] that shows the real tumour radius and the one obtained after 60 iterations. 
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Figure 13: Evolution of the tumour radius for each iteration: blue line corresponds to model-generated data without 
noise, and black line to model-generated data with a 5% of random noise. Red line represents the real observed 
radius at time T. 

First of all, Algorithm [7j] stops by a different reason: the functional can reach lower values in the case without 
noise than in the other one, so the derivative of the functional becomes an important stopping criteria. However, it 
is worth stressing that even with the presence of noise in the data, the method let us compute the parameters with 
a small error, in some cases of about 1%. Figure [13] shows the evolution of the tumour radius at time T for each 
iteration, comparing both cases with the real radius. 

9 Conclusions. 

The scientific community agrees that life's sciences, like biology or medicine, need the development of new tools 
in order to build models able to reproduce and to predict real phenomena. Over the last decades, a number of 
mathematical models for cancer onset and growth have been proposed QUI, and it became clear that these models 
are expected to success if the parameters involved in the modeling process are known. Or eventually, taking into 
account that some biological parameters may be unknown (especially in vivo), the model can be used to obtain 
them (2 ED. 

This paper, as already mentioned in Section 1, aims at offering a mathematical tool for the obtention of phe- 
nomenological parameters which can be identified by inverse estimation, by making suitable comparisons with 
experimental data. The inverse problem was stated as a PDE-constrained optimization problem, which was solved 
by using the adjoint method. The adjoint-based technique - although mathematically more complicated that the 
pattern search method used in lfl6l - has shown to work more efficiently, obtaining the results with better accuracy 
and with a less expensive numerical resolution. In addition, the gradient of the proposed functional is obtained and 
can be extended, in principle, to any number of unknown parameters. 

Focusing on further developments of the mathematical tools, it is worth mentioning that the numerical reso- 
lution proposed in this paper is in some aspects challenging and several numerical procedures were introduced in 
order to deal with non-linearities and singularities in the adjoint system of PDEs. 
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In addition, we remark that the parameter estimation via PDE-constrained optimization is a general approach 
that can be used, for instance, to consider the effects of chemotherapy. We are currently working in this line and 
also in the resolution of the optimization problem but after discretizing the original system of PDEs, in order to 
compare and contrast the performance of both methods. 
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Appendix 

Obtaining the adjoint problem. 



In this section we show the calculations involved in order to obtain the adjoint equations (|5.4| )- (|5.17| ). As stated in 
Section|4j the adjoint equations constitute a system of PDEs, with unknown X, given by (14.41) . To ease calculations, 
let us fix g = [n, v,c,s] T G t y. Hence, by ( 13.41 ), we have that 

a? =w Jo I l N M- N *M\ n M dtd y+& j o [s(t)-s*(t)] s (t)dt. gd 

On the other hand, (|f is obtained by using (15 - lb - In what follows, we shall obtain equivalent expressions 

for each of the nine terms of the summation (~jzg,X), which are associated with the nine constraints given by E in 
(1331 



Constraint 1 Consider the expression 

ff 

Jo Jo 



V-yS' s'S-S's at AT vS-Vs , . 
n t H 77— n y - — N y y + N y — ^ {a- bN) n 



S 2 -y s i 

- N (^ c -^ Nc ~ bn 



Xidtdy. 



Using integration by parts repeatedly and the facts that V(l,t) = S'(t) and V(0,t) = 0, it yields 



1 r T 



If 

Jo Jo 



S' S' V- V 

-ht + -jh + y^hy ~ y^l - -j^Xu, -(a- 2bN) Xi 

+ / / — Xivdtdy 
Jo Jo »j 



ndtdy 



f l f T (da db \, 

-jol N {^-^ N h cdtdy 

+ Jo Jo {-Y^ + ^-^ijsdtdy 
+ f (Xi (y, T)n(y, T) - X x (y,0)#i(y,0)) dy 

_ siT)i: mn^ii dy . 



(-2) 



Constraint 2 We have the following expression 

r 1 r T ( 2 db \ „ 

/ / v y + -v - ^NSc - bSn - bNs X 2 dtdy. 
Jo Jo \ y dC J 

In this case we have to integrate by parts just in the first term, because it is the only one that has a derivative of 
g, in this case v y . So we obtain that the second term in the inner product is 



o Jo 



2 \ db 
-X 2y + -X 2 I v — —NSX 2 c — bSX 2 n — bNX 2 s 
y J aC 



dtdy 



+ 



(-3) 



(X 2 (l,t)v(\,t)-X 2 (0,t)v(0,t))dt. 
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Constraint 3 The expression to be taken into account is 



lo fo i Cyy + \ Cy ~ kslyi ~ a% Ns2 ° ~ 2kNSs ^) ^ dtd y 



Because of the presence of second order derivatives we shall perform integration by parts twice. First of all, 
note that 

2 13 2 ^ 

Cyy + ~Cy = - — j (yc) . 

y y ay 1 

Then, we have, 



[ (cyy + lc^dy = f Q J^(yc)dy 



ho 

— a-OKO 

y dy 



y=l d /h\ v=1 r l d 2 (x 3 

-yc 



dy\y 



+ k yC 3f\l W 



To evaluate the limits, we assume that X 3 (y,t)c(y,t) — > when y — > and apply l'Hopital's rule. Hence, the 
third term is equal to 

fi [T ( 2X 3v 2X 3 \ ( , dk , \ 
J 1 i V ~ + ~^ J ° ~ \ + 3C + 2yUV5 * ) ^ 3 

+ £ (2X 3 (l,t)c(l,t) + X 3 (l,t)c y (l,t) -X 3y (l,t)c(l,t)) dt (.4) 



T 



(3h(0,t)c y (0,t) + hy(0,t)c(0,t)) dt. 



Constraint 4 The corresponding term in the inner product is 

J 1 (v(\,t)-s'(t))U(t)dt, 
where there is just one derivative of g involved. Integrating by parts we obtain: 



f (X 4 {t)v(l,t) + l 4t (t)s(t))dt-X 4 {T)s(T) + X 4 {0)s{0). (.5) 
Jo 



/o 

Constraint 5 Because this term is free of derivatives, there is nothing to do with it, remaining: 

f 

Jo 



X 5 (t)v(0,t)dt. (.6) 



Constraint 6 In this case, again, the corresponding term in the inner product remains: 







X 6 (t)c(l,t)dt. (.7) 



Constraint 7 Even though this term has a derivative, c y , it remains unchanged because function depends only 
on time: 



[ 7 \ 7 (t)cJ0,t)dt. (.8) 
Jo 
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Constraint 8 Once more, because of the lack of derivatives the term remains unchanged: 

[\ 8 (y)n(y,0)dy. (.9) 
Jo 

Constraint 9 In this case, the term is just the product of two real numbers: 

V(o). (.10) 

Obtaining the adjoint equations Note that according to equation (14.4b . we have to find X such that 

dJ fdE \\ 

where this equation should be valid for any direction g € < y. 

So, putting together equations ([j} with (T2l)- (1.10l) and choosing the directions conveniently, we get the system 
of equations which constitutes the adjoint problem (15.4I )- (15.17I ). 
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